'''
Introduction
============

`pyearthquake` is a light-weight module with basic earthquake engineering utilities.
The aim of this module is to help engineers play with earthquakes by python programming.

*Features*
--------

1. Convert .AT2 file downloaded from PEER database to nicely aligned form.  
2. Select and match ground motions to the target spectrum with Monte Carlo simulation and Least Square optimization.
3. Beautify the ground motions by truncating, trimming, changing dt, and renaming.
4. Store a suite of motions in a single file or separate files.
5. Generate displacement spectrums and pesudo acceleration spectrums.
6. Inspect a suite or a motion with nice html interactive plots, thanks to [plotly](https://plot.ly).  
7. Create Code Spectrums. Currently support Chinese Code.

Fact Sheet  
----------  
Version: 1.0   
Author: [Hanlin Dong](http://www.hanlindong.com)  
License: MIT  
Last Update: 2020-04-08   

Content Briefing
----------------

Currently, it has three main classes: `Spectrum`, `Motion`, and `Suite`

`Spectrum` : a class storing seismic spectrum.  
    Some static methods are created confining to Code Spectrums.

`Motion` : a class storing acceleration timehistory.  
    Motion has utility functions to nicely parse PEER motion file.  
    Motion also generates spectrums.

`Suite` : a class storing a suite of motions and their target spectrum.  
    Suite has utilities to load a PEER file folder.  
    Suite has optimization functions to select a bunch of earthquakes out of a large number of earthquakes
        with Monte Carlo simulation and bounded least square optimization.
        The advantage is that the shape of the accels is not changed unlike wavelet methods,
        but the average match result is still good.  
    Suite can write all the acceleration files seperately or integrally.   
    Suite can plot all the accelerations and spectrums interactively for users to inspect.  
    Suite can generate .tcl file to be used in OpenSees.  

All the classes mentioned above can be stored as string file.
    Currently they are json files with ext: .spectrum, .motion, .suite.
    These files can be read by static methods in every class.

Change log:
-----------
2020-04-08 13:53:32 v1.0  

- Join the three classes to one single file, refactor some functions.  
- Unify some methods, and create good documentation.
- Add plotly interactive plots.
'''

import os
import re
import json
import random
import numpy as np
import matplotlib.pyplot as plt
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from scipy.signal import lsim
from scipy.interpolate import interp1d
from scipy.optimize import least_squares


class Spectrum():
    '''A class storing a spectrum.  
    It can be init-ed simply with (periods, spectrum) sampling points.  
    It can also be init-ed by static functions to generate Code-confined spectrums.

    Parameters  
    ----------  
    `periods` : List, numpy.ndarray  
        List of sampling periods.  
        It can be generated by the static function `Spectrum.periods_standard` .
        Otherwise, users can specify desired period sampling points.  

    `spectrum` : List, numpy.ndarray  
        List of the spectrum points corresponding to sampling periods.  
    '''

    def __init__(self, periods, spectrum):
        self.periods = np.array(periods)
        '''numpy.ndarray, sampling periods'''
        self.spectrum = np.array(spectrum)
        '''numpy.ndarray, spectrum corresponding to sampling periods.'''

    def match(self, target, start, end, step=0.1):
        '''Match the current spectrum with a target spectrum. Use Least Square method.  

        Parameters  
        ----------  
        `target` : `Spectrum`  
            the spectrum to be matched.  

        `start`, `end` : float  
            the matching period boundaries.  

        `step` : float
            evenly sample the periods with the step.

        Returns
        -------
        `coefficient` : float  
            the coefficient to match the current spectrum to the target.

        `srss` : float   
            the srss after match.
        '''
        periods = np.arange(start, end+step, step)
        this_spectrum = np.array([self.get(p) for p in periods])
        that_spectrum = np.array([target.get(p) for p in periods])
        coefficient = sum(this_spectrum * that_spectrum) / sum(this_spectrum ** 2)
        srss = (np.sum((this_spectrum * coefficient - that_spectrum) ** 2)) ** 0.5
        return coefficient, srss

    def get(self, period):
        '''Get the spectrum value. Either a list or a number is accepted.  
        If the period exceeds the largest sampling period, the last spectrum value is returned.  
        The spectrum value not on sampling points are linearly interpolated.  

        Parameters
        ----------
        `period` : number, list, numpy.ndarray  
            The interested period.

        Returns
        -------
        `spectrum` : number, numpy.ndarray  
            if period is list-like, return a ndarray.
            if period is a number, return a number.
        '''
        if not hasattr(period, '__iter__'):
            i = self.periods.searchsorted(period)
            if self.periods[i] == period:
                return self.spectrum[i]
            elif i >= len(self.periods):
                print("Warning: given period exceeds the sampling range.")
                return self.spectrum[-1]
            else:
                return self.spectrum[i] - (self.spectrum[i]-self.spectrum[i-1]) / (self.periods[i] - self.periods[i-1]) * (self.periods[i] - period)
        else:
            return np.array([self.get(p) for p in period])

    def scale(self, factor_x=1, factor_y=1):
        '''Scale the spectrum.

        Parameters
        ----------
        `factor_x` : number  
            scaling factor in x direction (periods)  

        `factor_y` : number  
            scaling factor in y direction (spectrum)

        Returns
        -------
        `spectrum` : `Spectrum`
        '''
        return Spectrum(periods=self.periods*factor_x, spectrum=self.spectrum*factor_y)

    def plot(self, filename=None, engine='pyplot'):
        '''Plot the spectrum.

        Parameters
        ----------
        `filename` : string, None  
            if None, show it. Else, save it to the filename  

        `engine` : 'pyplot' or 'plotly'  
            Plot engine.

        Outputs
        -------
        A file with name `filename`.png or `filename`.html will be generated.
        '''
        if engine == 'pyplot':
            plt.figure(figsize=(4, 3))
            plt.plot(self.periods, self.spectrum, marker='o')
            plt.xlabel('Period (s)')
            plt.ylabel('Spetrum')
            plt.xlim([0, self.periods[-1]])
            if filename is None:
                plt.show()
            else:
                plt.savefig(filename+'.png', dpi=200)
        elif engine == 'plotly':
            fig = go.Figure()
            fig.add_scatter(dict(
                x='self.periods',
                y='self.spectrum',
                mode='lines+markers',
            ))
            if filename is None:
                fig.show()
            else:
                fig.write_html(filename+'.html')
        else:
            print('ERROR: engine should be in "pyplot" or "plotly"')

    def to_dict(self):
        '''Convert the instance properties to a dict.  

        Returns
        -------
        `it` : dict
        '''
        return dict(
            periods=self.periods.tolist(),
            spectrum=self.spectrum.tolist()
        )

    def to_json(self):
        '''Convert the properties to json.  

        Returns
        -------
        `json` : string  
        '''
        return json.dumps(self.to_dict())

    def save(self, filename):
        '''Save the spectrum to a file. Currently use json form.  

        Parameters
        ----------
        `filename` : string  
            The file name to be saved.   
            The extension '.spectrum' will be automatically added.  

        Outputs
        -------
        A file with name `filename`.spectrum will be generated.
        '''
        with open(filename + '.spectrum', 'w') as f:
            f.write(self.to_json())

    def write_csv(self, filename):
        '''Write the spectrum to a csv file.  

        Parameters
        ----------
        `filename` : string  
            The extension '.csv' will be automatically added.

        Outputs
        -------
        A file named `filename`.csv will be generated.
        '''
        np.savetxt(filename+'.csv', np.vstack(
            [self.periods, self.spectrum]).T, fmt="%.3f,%.7f", header='Period,Spectrum')

    @staticmethod
    def read_csv(filename):
        '''Create an instance by reading the written csv file.  

        Parameters
        ----------
        `filename` : string  
            The filename to read.

        Returns
        -------
        `spectrum` : `Spectrum`  
            New instance.
        '''
        if not filename.endswith('.csv'):
            filename += '.csv'
        data = np.loadtxt(filename, delimiter=',')
        return Spectrum(periods=data[:, 0], spectrum=data[:, 1])

    @staticmethod
    def read_dict(it):
        '''Create an instance from the dict.  

        Parameters
        ----------
        `it` : dict  
            The dict generated by `Spectrum.to_dict`

        Returns
        -------
        `spectrum` : `Spectrum`  
            New instance.
        '''
        return Spectrum(periods=it['periods'], spectrum=it['spectrum'])

    @staticmethod
    def read_json(string):
        '''Create an instance from json string.  

        Parameters
        ----------
        `string` : str  
            The json string 

        Returns
        -------
        `spectrum` : `Spectrum`  
            New instance.
        '''
        it = json.loads(string)
        return Spectrum.read_dict(it)

    @staticmethod
    def load(filename):
        '''Load the spectrum from file.  

        Parameters
        ----------
        `filename` : str  

        Returns
        -------
        `spectrum` : `Spectrum`  
            New instance.
        '''
        if not filename.endswith('.spectrum'):
            filename += '.spectrum'
        with open(filename, 'r') as f:
            string = f.read()
        return Spectrum.read_json(string)

    @staticmethod
    def periods_standard(until=6):
        '''Get a standard sampling of the periods.  
        The sampling points are more sparse with the increase of period.  

        Parameters
        ----------  
        `until` : number  
            The end of the period sampling list. It should >= 4.0.

        Returns
        -------
        `periods` : numpy.ndarray  
            The standard sampling points.
        '''
        assert until >= 4
        return np.hstack([np.arange(0, 0.2, 0.01),
                          np.arange(0.2, 0.5, 0.02),
                          np.arange(0.5, 1.0, 0.05),
                          np.arange(1.0, 2.0, 0.1),
                          np.arange(2.0, 4.0, 0.2),
                          np.arange(4.0, until+0.5, 0.5)])

    @staticmethod
    def chinese_point(period, intensity, major, site, group, damping=0.05):
        '''Get spectrum point for Chinese code, given a period.  

        Parameters
        ----------  
        `period` : number  
            The given period.  

        `intensity` : number  
            Intensity represented by a number.   
            Use 7.5 or 8.5 for half-level intensity.  

        `major` : boolean
            If True, it's for major (rare) earthquake.  
            If False, for minor (frequent) earthquake.  

        `site` : int
            Site number (Roman number in the code. 0 for I0, 1 for I1.)  

        `group` : int  
            Site class (Chinese number in the code.)  

        `damping` : float  
            Spectrum damping.  

        Returns
        -------  
        `value` : float
            the specific spectrum point value corresponding to the period.  
        '''
        # alpha max
        if intensity == 6:
            alphamax = 0.28 if major else 0.04
        elif intensity == 7:
            alphamax = 0.50 if major else 0.08
        elif intensity == 7.5:
            alphamax = 0.72 if major else 0.12
        elif intensity == 8:
            alphamax = 0.9 if major else 0.16
        elif intensity == 8.5:
            alphamax = 1.2 if major else 0.24
        elif intensity == 9:
            alphamax = 1.40 if major else 0.32
        else:
            raise("ERROR: intensity not in [6, 7, 7.5, 8, 8.5, 9]")
        # tg
        tg_table = [[0.20, 0.25, 0.35, 0.45, 0.65],
                    [0.25, 0.3, 0.4, 0.55, 0.75],
                    [0.3, 0.35, 0.45, 0.65, 0.9]]
        assert group in [1, 2, 3] and site in [0, 1, 2, 3, 4]
        tg = tg_table[group-1][site]
        # other parameters
        gamma = 0.9 + (0.05 - damping) / (0.3 + 6*damping)
        eta1 = 0.02 + (0.05 - damping) / (4 + 32*damping)
        eta2 = 1 + (0.05 - damping) / (0.08 + 1.6*damping)
        # calculate
        if period < 0.1:
            return ((10*eta2 - 4.5) * period + 0.45) * alphamax
        elif period < tg:
            return eta2 * alphamax
        elif period < 5*tg:
            return (tg/period)**gamma * eta2 * alphamax
        elif period <= 6:
            return (eta2 * 0.2**gamma - eta1 * (period-5*tg)) * alphamax
        else:
            return (eta2 * 0.2**gamma - eta1 * (6-5*tg)) * alphamax

    @staticmethod
    def chinese(intensity=8, major=True, site=3, group=2, damping=0.05):
        '''Return a new spectrum instance that confines to Chinese Code.  
        Use `Spectrum.periods_standard` to get sampling periods.  

        Parameters
        ----------
        refer to `Spectrum.chinese_point`

        Returns
        -------
        `spectrum` : `Spectrum`
        '''
        periods = Spectrum.periods_standard(until=6)
        spectrum = [Spectrum.chinese_point(p, 
                                           intensity=intensity,
                                           major=major,
                                           site=site,
                                           group=group,
                                           damping=damping) for p in periods]
        return Spectrum(periods, spectrum)


class Motion(object):
    '''A class plays with ground motion. Store accel time history and accel spectrum only.  
    '''

    def __init__(self, folder, name, damping=0.05):
        '''
        Parameters
        ----------
        `folder` : str   
            the folder this motion is stored.

        `name` : str  
            The motion name. The motion will be saved to '{folder}/{name}.motion'

        `damping` : float  
            The damping ratio of the motion spectrum.  
        '''
        self.folder = folder
        '''str. See `Motion`.'''
        os.makedirs(folder, exist_ok=True)
        self.name = name
        '''str. See `Motion`.'''
        self.accel = None
        '''numpy.ndarray. Acceleration time history with evenly distributed sampling points.'''
        self.dt = None
        '''float. Acceleration time sampling delta t.'''
        self.information = ''
        '''str. Store massive infomation about the motion. '''
        self.spectrum = None
        '''`Spectrum`. The spectrum of the motion.'''
        self.damping = damping
        '''float. The damping ratio of the spectrum. Note: init with 0.05.'''

    def get_PGA(self):
        '''
        Returns
        -------
        `value` : float  
            The peak ground acceleration.  
        '''
        return np.max(abs(self.accel))

    def get_duration(self):
        '''
        Returns
        -------
        `value` : float  
            The duration of the motion.  
        '''
        return (len(self.accel) - 1) * self.dt

    def get_times(self):
        '''
        Returns
        -------
        `value` : numpy.ndarray  
            The time series of the motion  
        '''
        return np.arange(len(self.accel)) * self.dt

    def generate_spectrum(self, periods=None):
        '''Use scipy LTS signal processing module to generate the spectrum.  
        Note that the accel spectrum is pseudo-spectrum derived from disp spectrum.  
        So it's only accurate when damping is small.  
        If `Motion.dt` is still None, which usually occurs if the PEER file is not well formated, 
            a blank spectrum will be returned.  

        Parameters
        ----------
        `periods` : List-like, None    
            The period sampling points.   
            If is None, `Spectrum.periods_standard` will be called.

        Returns
        -------
        `spectrum` : `Spectrum`  
            The generated spectrum.

        Outputs
        -------
        `Motion.spectrum` will be changed.
        '''
        if periods is None:
            periods = Spectrum.periods_standard()
        if self.dt is not None:
            print("Generating spectrum ...")
            omegas = 2 * np.pi / periods[1:]
            num = np.array([-1])
            spectrum = np.zeros(len(periods))
            spectrum[0] = self.get_PGA()
            for i, omega in enumerate(omegas):
                den = np.array([1, 2*self.damping*omega, omega**2])
                _, y, _ = lsim((num, den), self.accel, self.get_times())
                spectrum[i+1] = np.max(abs(y)) * omega ** 2
            self.spectrum = Spectrum(periods, spectrum)
        else:
            print("Dt is not known. Assigning empty spectrum.")
            self.spectrum = Spectrum(periods=[], spectrum=[])
        return self.spectrum

    def load_peer(self, peer_filename):
        '''Load acceleration time-history files downloaded from PEER.  
        The strings in the file will be parsed and also stored in self.information.  

        Parameters
        ----------
        `peer_filename` : str

        Outputs
        -------
        `Motion.dt`, `Motion.accel`, `Motion.information`, `Motion.spectrum` will all be changed.  
        '''
        self.information += f'peer_filename: {peer_filename}'
        try:
            f = open(peer_filename, 'r')
            for _ in range(3):
                self.information += f.readline()
            s = f.readline()
            self.information += s
            re_dt = re.compile(r".*[ 0](\.[0-9]*) ?SEC")
            mt = re_dt.match(s)
            try:
                self.dt = float(mt.group(1).replace('.', '0.'))
                print(f'Extracted dt = {self.dt:.4f}')
            except:
                print(f'ERROR: cannot find dt from {s}')
                self.dt = 0.000001
            series = []
            datastring = f.readlines()
            for dtline in datastring:
                for dt in dtline.split("  "):
                    if dt:
                        try:
                            series.append(float(dt.replace(".", "0.")))
                        except ValueError:
                            break
            self.accel = np.array(series)
            print(f'Extracted npts = {len(self.accel)}')
        except IOError:
            print(f'No file named {self.name} is found.')

    def convert_peer(self, peer_filename):
        '''A recipe to convert the peer acceleration file to the .motion file.  

        Parameters
        ----------
        `peer_filename` : str  

        `folder` : str  
            the folder to save .motion file.

        Outputs
        -------
        A .motion file will be created in `folder` with the PEER filename.
        '''
        print(f'Converting PEER file {peer_filename} to folder {self.folder} with name {self.name}')
        self.load_peer(peer_filename)
        self.generate_spectrum()
        self.save()

    def trim(self, new_name, new_folder=None, left=1.0e-3, right=1.0e-3, prepend_zero=True, log=True, log_folder=None):
        '''Trim the motion with two coefficients of PGA.  
        This method is used when a long near-zero time-history is stored in the data file.  

        Parameters
        ----------
        `new_name` : str  

        `new_folder` : str, None  

        `left`, `right` : float 
            The coefficient of PGA to trim the time history.
            Search from left or right, when the target point is found,
            search for the nearest near-zero point. Then trim off the small values.  


        `prepend_zero` : boolean  
            If True, prepend a zero value to the trimmed acceleration.

        `log` : boolean  
            If the log picture will be generated.

        `log_folder` : str, None
            Path for the folder to save the trim logs.
            If None, save in the previous folder.

        Returns
        -------
        `motion` : `Motion`  
            The trimmed motion.

        Outputs
        -------
        If log, plot a figure in log_folder named '{new_name}_trimlog.png'
        '''
        mo = self.copy()
        mo.name = new_name
        mo.folder = new_folder if new_folder is not None else self.folder
        pga = self.get_PGA()
        # search left
        i = 0
        while abs(self.accel[i]) < pga * left:
            i += 1
        while i > 0 and self.accel[i] * self.accel[i-1] > 0:
            i -= 1
        i_left = i
        # search right
        i = -1
        while abs(self.accel[i]) < pga * right:
            i -= 1
        while i < -1 and self.accel[i] * self.accel[i+1] > 0:
            i += 1
        i_right = len(self.accel) + i
        if prepend_zero:
            mo.accel = np.hstack([0, self.accel[i_left:i_right]])
        else:
            mo.accel = self.accel[i_left:i_right]
        mo.generate_spectrum()
        print(
            f'Motion is trimmed from {self.name}:{len(self.accel)*self.dt:.3f}s to {mo.name}:{len(mo.accel)*mo.dt:.3f}s.')
        if log:
            log_folder = self.folder if log_folder is None else log_folder
            os.makedirs(log_folder, exist_ok=True)
            _, ax = plt.subplots(3, 1, figsize=(7, 5))
            ax[0].plot(self.get_times(), self.accel, label='before')
            ax[1].plot(mo.get_times(), mo.accel, label='after')
            ax[2].plot(self.spectrum.periods, self.spectrum.spectrum, label='before')
            ax[2].plot(mo.spectrum.periods, mo.spectrum.spectrum, label='after')
            for i in range(3):
                ax[i].legend()
            plt.suptitle(f'Trim log: {self.name}:{len(self.accel)*self.dt:.3f}s to {mo.name}:{len(mo.accel)*mo.dt:.3f}s.')
            plt.savefig(f'{log_folder}/{mo.name}_trimlog.png', dpi=200)
        return mo

    def change_dt(self, dt, new_name, new_folder=None, log=True, log_folder=None):
        '''Change dt for the accel. Linearly interpolate the acceleration, and re-sample.
        Note: no low-pass filters are used. May cause new peaks in the Fourier spectrum.  

        Parameters
        ----------
        `dt` : float  
            New sampling delta t. 

        `new_name` : str  

        `new_folder` : str, None  

        `log` : boolean  
            If log, create a log picture comparing the before and after accels and spectrum.

        `log_folder` : str, None  
            Path for the folder to save the trim logs.
            If None, save in the previous folder.

        Returns
        -------
        `motion` : `Motion`  

        Outputs
        -------
        If log, create a picture named '{log_folder}/{new_name}_dtlog.png'.
        '''
        if dt == self.dt:
            print(f'Motion {self.name} do not need to change dt')
            return self.copy()
        else:
            mo = self.copy()
            mo.name = new_name
            mo.folder = new_folder if new_folder is not None else self.folder
            f = interp1d(self.get_times(), self.accel)
            npts = np.floor(self.dt * len(self.accel) / dt)
            mo.accel = f(np.arange(npts)*dt)
            mo.dt = dt
            mo.generate_spectrum()
            print(f'Motion dt changed from {self.name}:{self.dt:.4f}s to {mo.name}:{dt:.4f}s.')
        if log:
            log_folder = self.folder if log_folder is None else log_folder
            os.makedirs(log_folder, exist_ok=True)
            _, ax = plt.subplots(3, 1, figsize=(7, 5))
            ax[0].plot(self.get_times(), self.accel, label='before')
            ax[1].plot(mo.get_times(), mo.accel, label='after')
            ax[2].plot(self.spectrum.periods, self.spectrum.spectrum, label='before')
            ax[2].plot(mo.spectrum.periods, mo.spectrum.spectrum, label='after')
            for i in range(3):
                ax[i].legend()
            plt.suptitle(f'Change dt: {self.name}:{self.dt:.4f}s to {mo.name}:{dt:.4f}s.')
            plt.savefig(f'{log_folder}/{mo.name}_dtlog.png', dpi=200)
        return mo

    def truncate(self, time, new_name, new_folder=None, log=True, log_folder=None):
        '''Truncate the acceleration. 
        Sometimes, a ground motion file may have two peaks or more. 
        So if only one peak is needed, this method can be used to skip the other peaks.  

        Parameters
        ----------
        `time` : float
            The accels after the time will be cut off.  

        `new_name` : str

        `new_folder` : str, None  

        `log` : boolean  

        `log_folder` : str, None  

        Returns
        -------
        `motion` : `Motion`

        Outputs
        -------
        If log is True, a figure named '{log_folder}/{new_name}_trunclog.png' will be created.
        '''
        i = round(time / self.dt)
        mo = self.copy()
        mo.name = new_name
        mo.folder = new_folder if new_folder is not None else self.folder
        if i >= len(self.accel):
            print(f"WARNING: time is longer than duration. {self.name} is not truncated.")
            return mo
        mo.accel = self.accel[:i]
        mo.generate_spectrum()
        print(f'Motion truncated: {self.name}:{len(self.accel)*self.dt:.3f}s to {mo.name}:{len(mo.accel)*mo.dt:.3f}s.')
        if log:
            log_folder = self.folder if log_folder is None else log_folder
            os.makedirs(log_folder, exist_ok=True)
            _, ax = plt.subplots(3, 1, figsize=(7, 5))
            ax[0].plot(self.get_times(), self.accel, label='before')
            ax[1].plot(mo.get_times(), mo.accel, label='after')
            ax[2].plot(self.spectrum.periods, self.spectrum.spectrum, label='before')
            ax[2].plot(mo.spectrum.periods, mo.spectrum.spectrum, label='after')
            for i in range(3):
                ax[i].legend()
            plt.suptitle(f'Truncate: {self.name}:{len(self.accel)*self.dt:.3f}s to {mo.name}:{len(mo.accel)*mo.dt:.3f}s.')
            plt.savefig(f'{log_folder}/{mo.name}_trunclog.png', dpi=200)
        return mo

    def scale(self, new_name, new_folder=None, factor_x=1, factor_y=1):
        '''Scale the ground motion.  

        Parameters
        ----------
        `new_name` : str

        `new_folder` : str, None

        `factor_x` : float  
            The scaling factor on the time/period axis.  

        `factor_y` : float
            The scaling factor on the accel axis.

        Returns
        -------
        `motion` : Motion
        '''
        mo = self.copy()
        mo.name = new_name
        mo.folder = new_folder if new_folder is not None else self.folder
        if factor_x != 1:
            mo.dt = self.dt * factor_x
        if factor_y != 1:
            mo.accel = self.accel * factor_y
        mo.spectrum = self.spectrum.scale(factor_x, factor_y)
        return mo

    def plot_timehistory(self, folder=None, save=True, engine='pyplot'):
        '''Plot the time history of the acceleration.

        Parameters
        ----------
        `folder` : str, None  
            Store folder. If None, use `Motion.folder` 

        `save` : boolean
            Save the figure or not. 

        `engine` : 'pyplot' or 'plotly'  

        Outputs
        -------
        if save is True, create a figure '{folder}/{self.name}_th'.
        '''
        folder = self.folder if folder is None else folder
        if engine == 'pyplot':
            plt.figure()
            t = self.get_times()
            plt.plot(t, self.accel)
            if save:
                plt.savefig(f'{folder}/{self.name}_th.png', dpi=200)
            else:
                plt.show()
        elif engine == 'plotly':
            fig = go.Figure()
            fig.add_trace(go.Scatter(
                x=self.get_times(),
                y=self.accel,
                mode='lines',
                name=self.name
            ))
            if save:
                fig.write_html(f'{folder}/{self.name}_th.html')
            else:
                fig.show()
        else:
            print('ERROR: engine should be in "pyplot" or "plotly"')

    def plot(self, folder=None, save=True, engine='pyplot'):
        '''Plot the time history and the spectrum.

        Parameters
        ----------
        See `Motion.plot_timehistory`

        Outputs
        -------
        Create a figure named '{folder}/{name}'
        '''
        folder = self.folder if folder is None else folder
        if engine == 'pyplot':
            _, ax = plt.subplots(2, 1, figsize=(7, 8))
            ax[0].plot(self.get_times(), self.accel)
            ax[0].set_xlabel('Time(s)')
            ax[0].set_ylabel('Acceleration(g)')
            ax[0].set_title(f'PGA={self.get_PGA():.4f}')
            ax[1].plot(self.spectrum.periods, self.spectrum.spectrum)
            ax[1].set_xlabel('Period(s)')
            ax[1].set_ylabel('Spectrum Acceleration(g)')
            ax[1].set_title(f'damping={self.damping:.4f}')
            plt.subplots_adjust(hspace=0.3)
            if save:
                plt.savefig(f'{folder}/{self.name}.png', dpi=200)
            else:
                plt.show()
        elif engine == 'plotly':
            fig = make_subplots(rows=1, cols=2)
            fig.add_trace(go.Scatter(
                x=self.get_times(),
                y=self.accel,
                mode='lines',
                name='accel'
            ), row=1, col=2)
            fig.add_trace(go.Scatter(
                x=self.spectrum.periods,
                y=self.spectrum.spectrum,
                mode='lines',
                name='spectrum'
            ), row=1, col=1)
            if save:
                fig.write_html(f'{folder}/{self.name}.html')
            else:
                fig.show()
        else:
            print('Engine should be in "pyplot" or "plotly"')

    def write_accel(self, filename, time=True):
        '''Write the acceleration timehistory to a single file.  

        Parameters
        ----------
        `filename` : str  
            The file path to store the accelration.  

        `time` : boolean  
            If True, a time column is created.
            To know dt in the file content, time should be True.
            Otherwise, dt should be stored in the filename.

        Outputs
        -------
        An text file with acceleration information will be written.
        '''
        if time:
            data = np.vstack([self.get_times(), self.accel])
            np.savetxt(filename, data.T, fmt="%.4f,%.7e", header="Time,Accel")
        else:
            np.savetxt(filename, self.accel, fmt="%.7e")

    def write_info(self, filename):
        '''Write the information to a single file.  

        Parameters
        ----------
        `filename` : str

        Outputs
        -------
        A text file with information will be created.
        '''
        with open(filename, 'w') as f:
            f.write(f'dt={self.dt:.3}\n')
            f.write(f'npts={len(self.accel)}\n')
            f.write(f'damping={self.damping:.3}\n')
            f.write(self.information)

    def write_spectrum(self, filename):
        '''Write the spectrum to a single file. Use the form of csv.

        Parameters
        ----------
        `filename` : str

        Outputs
        -------
        A csv file will be created.
        '''
        self.spectrum.write_csv(filename)

    def read_accel(self, filename, dt=True):
        '''Load an single-column acceleration file to the current instance.

        Parameters
        ----------
        `dt` : float, boolean
            if dt is float, `Motion.dt` will be set to dt.
            Else if dt is True, `Motion.dt` will be found in the file.

        Outputs
        -------
        Change `Motion.dt` and `Motion.accel` in the CURRENT instance.
        '''
        if not isinstance(dt, float):
            data = np.loadtxt(filename)
            self.accel = data[:, 1]
            self.dt = data[1, 0]
        else:
            self.accel = np.loadtxt(filename)
            self.dt = dt

    def read_spectrum(self, filename):
        '''Read the spectrum from csv file.  

        Parameters
        ----------
        `filename` : str

        Outputs
        -------
        The `Motion.spectrum` of the CURRENT instance will be changed.
        '''
        spectrum = Spectrum.read_csv(filename)
        self.spectrum = spectrum

    def read_info(self, filename):
        '''Read info from text file.  

        Parameters
        ----------
        `filename` : str

        Outputs
        -------
        The `Motion.information` of the CURRENT instance will be changed.
        '''
        with open(filename, 'r') as f:
            self.dt = float(f.readline().replace("dt=", ""))
            _ = int(f.readline().replace("npts=", ""))
            self.damping = float(f.readline().replace("damping=", ""))
            self.information = f.read()

    def to_dict(self):
        '''Convert the instance properties to dict.  

        Returns
        -------
        `it` : dict
        '''
        return dict(
            dt=self.dt,
            npts=len(self.accel),
            damping=self.damping,
            information=self.information,
            accel=self.accel.tolist(),
            spectrum=self.spectrum.to_dict(),
            folder=self.folder,
            name=self.name,
        )

    def to_json(self):
        '''Convert the instance properties to string using json format.  

        Returns
        -------
        `string` : str
        '''
        return json.dumps(self.to_dict())

    def save(self):
        '''Save the instance to a file, which is defined by `Motion.folder` and `Motion.name`  

        Outputs
        -------
        A file named '{self.folder}/{self.name}.motion' will be created.
        '''
        with open(f'{self.folder}/{self.name}.motion', 'w') as f:
            f.write(self.to_json())

    @staticmethod
    def read_dict(it):
        '''Create a instance by the dict generated by `Motion.to_dict`.  

        Parameters
        ----------
        `it` : dict

        Returns
        -------
        `mo` : `Motion`
        '''
        mo = Motion(folder=it['folder'], name=it['name'])
        mo.dt = it['dt']
        mo.damping = it['damping']
        mo.information = it['information']
        mo.accel = np.array(it['accel'])
        mo.spectrum = Spectrum.read_dict(it['spectrum'])
        return mo

    @staticmethod
    def read_json(string):
        '''Create a new instance by the string generated by `Motion.to_json`.  

        Parameters
        ----------
        `string` : str

        Returns
        -------
        `mo` : `Motion`
        '''
        it = json.loads(string)
        return Motion.read_dict(it)

    @staticmethod
    def load(filename):
        '''Load the file created from `Motion.save` method.  

        Parameters
        ----------
        `filename` : str

        Returns
        -------
        `mo` : `Motion`
        '''
        if not filename.endswith('.motion'):
            filename += '.motion'
        with open(filename, 'r') as f:
            return Motion.read_json(f.read())

    def copy(self):
        '''Copy the instance.

        Returns
        -------
        `mo` : `Motion`
        '''
        return Motion.read_dict(self.to_dict())


class Suite(object):
    '''A class that contains mainly a series of `Motion` instances, 
    and corresponding series of factors, matching the suite of motions to the target `Spectrum`.  
    '''

    def __init__(self, folder, name, target_spectrum, period_bound):
        '''
        Parameters
        ----------
        `folder` : str  
            A path-like string for the suite.

        `name` : str
            A specific name for the suite. The instance will be saved in '{folder}/{name}.suite'  

        `target_spectrum` : `Spectrum`  
            The target spectrum of the suite.  

        `period_bound` : List<float> len=2  
            The matching period boundaries.
        '''
        self.folder = folder
        '''str. See `Suite`'''
        os.makedirs(folder, exist_ok=True)
        self.name = name
        '''str. See `Suite`'''
        self.motions = []
        '''List<`Motion`>. Store all the motions.'''
        self.factors = []
        '''List<float>. Store all the factors that should be applied to the motions.'''
        self.target_spectrum = target_spectrum
        '''`Spectrum`. See `Suite`'''
        self.period_bound = period_bound
        '''List<float>. See `Suite`'''
        self.period_samples = np.linspace(period_bound[0], period_bound[1], 200)
        '''numpy.ndarray. The period within the `period_bound` is evenly sampled, with 200 points.'''
        self.spectrum_matrix = None
        '''numpy.ndarray. A matrix with shape (len(`Suite.period_samples`), len(`Suite.motions`)).
        This matrix times the `Suite.factors` vector products the average spectrum with `Suite.period_samples`.'''
        self.target_vector = target_spectrum.get(self.period_samples)
        '''numpy.ndarray. The vector generated by the target_spectrum sampling the `Suite.period_samples`'''

    def set_folder(self, folder):
        '''Reset the folder for the suite and the motions.

        Parameters
        ----------
        `folder` : str

        Outputs
        -------
        `Suite.folder` and `Suite.motions` will be changed.
        '''
        if folder is not None:
            self.folder = folder
            os.makedirs(folder, exist_ok=True)
            for mo in self.motions:
                mo.folder = folder

    def add_motion(self, motion, factor=1):
        '''Add motion to the suite together with a factor, to keep the lengths are identical.

        Parameters
        ----------
        `motion` : `Motion`

        `factor` : number

        Outputs
        -------
        `Suite.motions` and `Suite.factors` will be changed.
        '''
        self.motions.append(motion)
        self.factors.append(factor)

    def load_motions_from_names(self, names, folder=None):
        '''Give a list of name and the folder, load motions.

        Parameters
        ----------
        `names` : List<str>  
            The list of motion names.

        `folder` : str, None  
            The folder of the motions.  
            If is None, use `Suite.folder`

        Outputs
        -------
        `Suite.motions` and `Suite.factors` will be changed.
        '''
        folder = self.folder if folder is None else folder
        for name in names:
            mo = Motion.load(f'{folder}/{name}.motion')
            self.add_motion(mo)
            print(f'Loaded motion {mo.name}')
        self.reset_spectrum_matrix()

    def load_motions_from_folder(self, folder=None):
        '''Find all the .motion file from the given folder, and load them.  

        Parameters
        ----------
        `folder` : str, None   
            The folder to search the .motion files.  
            If is None, use `Suite.folder` instead.  

        Outputs
        -------
        `Suite.motions` and `Suite.factors` will be changed.
        '''
        folder = self.folder if folder is None else folder
        filenames = os.listdir(self.folder)
        for filename in filenames:
            if filename.endswith('.motion'):
                mo = Motion.load(f'{folder}/{filename}')
                self.add_motion(mo)
                print(f'Loadad motion {mo.name}')
        self.reset_spectrum_matrix()

    def filter_by_IDs(self, ids, new_name, new_folder=None):
        '''Give a series of selected IDs, filter the motions.  

        Parameters
        ----------
        `ids` : List<int>  
            The list of indexes of the selected motions in the current suite.  

        `new_name` : str  

        `new_folder` : str, None

        Returns
        -------
        `suite` : `Suite`  
            A new suite with the selected motions and their factors only.  
        '''
        suite = self.copy()
        suite.name = new_name
        suite.set_folder(new_folder)
        suite.motions = [self.motions[i] for i in ids]
        suite.factors = [self.factors[i] for i in ids]
        suite.reset_spectrum_matrix()
        print(f'The motions have been filtered. {len(ids)} motions remain.')
        return suite

    def filter_by_file(self, filename, new_name, new_folder=None):
        '''Give a file with indexes selected.  

        Parameters
        ----------
        Refer to `Suite.filter_by_IDs`.  

        Returns
        -------
        Refer to `Suite.filter_by_IDs`.  
        '''
        with open(filename, 'r') as f:
            lines = f.readlines()
        ids = [int(l) for l in lines]
        return self.filter_by_IDs(ids=ids, new_name=new_name, new_folder=new_folder)

    def get_average_spectrum(self):
        '''Return a spectrum which is the average spectrum of the suite.

        Outputs
        -------
        `spectrum` : `Spectrum`  
        '''
        periods = self.motions[0].spectrum.periods
        average_spectrum = self.motions[0].spectrum.spectrum * self.factors[0]
        for i in range(1, len(self.motions)):
            average_spectrum += self.motions[i].spectrum.get(periods) * self.factors[i]
        average_spectrum /= len(self.motions)
        return Spectrum(periods, average_spectrum)

    def get_srss(self):
        '''Get the srss of the current averate spectrum to the target spectrum in the period boundary.

        Returns
        -------
        `srss` : float
        '''
        return np.linalg.norm(np.matmul(self.spectrum_matrix, self.factors) - self.target_vector).item()

    def reset_spectrum_matrix(self):
        '''Reset the spectrum matrix, which produces the average spectrum 
        in the boundary by multiplying factors vector.
        This method should be called everytime `Suite.motions` change.

        Outputs
        -------
        `Suite.spectrum_matrix` and `Suite.target_vector` will be changed.  
        '''
        mat = np.zeros((len(self.period_samples), len(self.motions)))
        for i in range(len(self.motions)):
            spectrum = self.motions[i].spectrum.get(self.period_samples)
            mat[:, i] = spectrum / len(self.motions)
        self.spectrum_matrix = mat
        self.target_vector = self.target_spectrum.get(self.period_samples)

    def eliminate_neg(self):
        '''Eliminate negative target vector by multiplying a factor for every individual.
        In other words, make the average spectrum above the target spectrum.

        Returns
        -------
        `factor` : float  
            The global factor that should be applied to each of the factors.  

        Outputs
        -------
        `Suite.factors` are all amplified by the returned `factor`.
        '''
        error = np.matmul(self.spectrum_matrix, self.factors) - self.target_vector
        argmin = np.argmin(error)
        factor = 1 / (1 + error[argmin] / self.target_vector[argmin])
        print(f'Multiply {factor:.4} to eliminate negative spectrum.')
        self.factors = [self.factors[i] * factor for i in range(len(self.factors))]
        return factor

    def match_individual_LSQ(self):
        '''Use Least square method for each individual motion to match the spectrum.

        Returns
        -------
        `residuals` : numpy.ndarray  
            srss of each motion matched to the target spectrum.  

        Outputs
        -------
        `Suite.factors` are changed.
        '''
        residuals = np.zeros(len(self.motions))
        for i in range(len(self.motions)):
            this = self.spectrum_matrix[:, i] * len(self.motions)
            that = self.target_vector
            factor = np.dot(this, that) / np.dot(this, this)
            residual = np.linalg.norm(this * factor - that)
            self.factors[i] = factor
            residuals[i] = residual
        print(f"Individual LSQ: Factor Max={np.max(self.factors):.4}, Residual Max={np.max(residuals):.4}")
        return residuals

    def filter_optimize(self, count, times, use_lsq=True, output_count=0, new_folder=None, lower_bound=0.6, upper_bound=1.4, dimension=100000):
        '''Use Monte Carlo simulation and Least Square method to optimize the suite.
        after running, a bunch of suites will be written to the new_folder.

        Parameters  
        ----------
        `count` : int  
            The number of motions in the resulting suite.  

        `times` : int  
            The number of times that MC simulation with dimension is run.  
            The total number of samples of MC simulation will be {times} * {dimension}.  

        `use_lsq` : boolean  
            Whether use LSQ optimizer to better optimize the suite.  

        `output_count` : int  
            The number of outputs. if 0, output all.  

        `new_folder` : str 
        `lower_bound`, `upper_bound` : float  
            For the purpose that the individual motions are not too far from the individually matched factors,
            the lower and upper boundaries of amplifying the individually matched factors are set.  
            This is only useful when `use_lsq` is set, to add boundaries to the LSQ optimizer.  

        `dimension` : int  
            Each time the Monte Carlo simulation is run, 
            means that {dimension} times of simulation have been done by using a matrix target.

        Returns
        -------
        List<`Suite`> with len=`output_count`.  
        The list is sorted. The best answer is the first result. However, visual check should be done.

        Outputs
        -------
        Every running after Monte Carlo simulation and LSQ optimization, the spectrum of the suite will be plotted.
        '''
        new_folder = new_folder if new_folder is not None else self.folder
        os.makedirs(new_folder, exist_ok=True)
        suites = []
        srsses = []
        for i in range(times):
            suite_new = self.filter_montecarlo(count=count, new_name=f'mc{i:03d}', new_folder=new_folder)
            suites.append(suite_new)
            srsses.append(suite_new.get_srss())
            if use_lsq:
                suite_new2 = suite_new.optimize_LSQ(
                    new_name=f'lsq{i:03d}', 
                    new_folder=new_folder, 
                    lower_bound=lower_bound, 
                    upper_bound=upper_bound
                )
                suites.append(suite_new2)
                srsses.append(suite_new2.get_srss())
        argsorted = np.argsort(np.array(srsses))
        results = [suites[i] for i in argsorted[:output_count]]
        for i, suite in enumerate(results):
            suite.name = f'no{i}-{suite.name}'
            suite.save()
            suite.plot_all_spectrums(engine='pyplot')
        return results

    def filter_montecarlo(self, count, new_name, new_folder=None, dimension=100000):
        '''Use Monte Carlo simulation to filter the suite to count numbers of motions.

        Parameters
        ---------- 
        Refer to `Suite.filter_optimize`.  

        Returns
        -------
        `suite` : `Suite`   
            The optimized result after {dimension} times of MC simulations.
        '''
        print(f'Running Monte Carlo simulation, dimension={dimension}')
        self.match_individual_LSQ()
        factors = np.array(self.factors)
        xs = np.zeros((len(self.motions), dimension))
        for i in range(dimension):
            indexes = random.sample(range(len(self.motions)), count)
            xs[indexes, i] = factors[indexes]
        mul = np.matmul(self.spectrum_matrix * len(self.motions) / count, xs)
        loss = mul - self.target_vector.reshape(-1, 1)
        argmin = np.argmin(loss, axis=0)
        lossmin = np.min(loss, axis=0)
        target_at_lossmin = self.target_vector[argmin]
        amps = 1 / (lossmin / target_at_lossmin + 1)
        lossAmp = mul * amps - self.target_vector.reshape(-1, 1)
        srss = np.linalg.norm(lossAmp, axis=0)
        arg = np.argmin(srss)
        factors = xs[:, arg] * amps[arg]
        suvives = np.where(factors > 1.0e-3)
        suite = self.copy()
        suite.name = new_name
        suite.set_folder(new_folder)
        suite.motions = [self.motions[i] for i in suvives[0]]
        suite.factors = [factors[i] for i in suvives[0]]
        suite.reset_spectrum_matrix()
        return suite

    def _loss(self, x):
        '''Utility function: calculate the loss for optimization.'''
        factors = x * np.array(self.factors)
        mul = np.matmul(self.spectrum_matrix, factors)
        loss = mul - self.target_vector
        argmin = np.argmin(loss)
        amp = 1 / (loss[argmin] / self.target_vector[argmin] + 1)
        return mul * amp - self.target_vector

    def optimize_LSQ(self, new_name, new_folder=None, lower_bound=0.6, upper_bound=1.4):
        '''Use Least Square method from scipy to optimize the factors.

        Parameters
        ----------
        Refer to `Suite.filter_optimize`.   

        Returns
        -------
        `suite` : `Suite`  
            Optimized suite.
        '''
        print('Running Bounded Least Square optimization ...')
        res = least_squares(self._loss,
                            np.ones(len(self.motions)),
                            bounds=(lower_bound, upper_bound),
                            verbose=0)
        factors = res.x * self.factors
        suite = self.copy()
        suite.name = new_name
        suite.set_folder(new_folder)
        suite.factors = factors.tolist()
        suite.eliminate_neg()
        return suite

    def load_peer_folder(self, peer_folder, h1=True, h2=True, v=False):
        '''From a folder load all the PEER motions into the suite.  
        Extract the downloaded peer zip file first.

        Parameters
        ----------
        `peer_folder` : str  
            The folder path of the extracted PEER data. Where `_SearchResults.csv` can be found.

        `h1`, `h2`, `v` : boolean  
            Whether load the three directional files or not.  

        Outputs
        -------
        `Suite.motions`, `Suite.factors` will be changed.  
        All the PEER motions will be saved in .motion format in the new folder.  
        '''
        with open(f'{peer_folder}/_SearchResults.csv', "r") as f:
            for _ in range(34):
                f.readline()
            i = 0
            positions = []
            if h1:
                positions.append(19)
            if h2:
                positions.append(20)
            if v:
                positions.append(21)
            line = f.readline()
            while line != '\n':
                cells = line.split(',')
                for pos in positions:
                    filename = cells[pos].replace(' ', '')
                    motion_name = filename.replace('.AT2', '')
                    print(f'Loading {motion_name}')
                    mo = Motion(folder=self.folder, name=motion_name)
                    mo.load_peer(peer_filename=f'{peer_folder}/{filename}')
                    mo.generate_spectrum()
                    mo.save()
                    self.add_motion(mo)
                    i += 1
                line = f.readline()
        self.reset_spectrum_matrix()
        print(f"Successfully loaded {i} earthquakes from peer folder {peer_folder} to {self.folder}.")
        return True

    def beautify(self, new_name, run_steps=[1,2,3,4], dt=0.01, new_folder=None, trunc_dict={}, left=1.0e-3, right=1.0e-3, prepend_zero=True):
        '''Beautify the suite in 4 steps:  

        1. truncate the motions, see `Motion.truncate`  
        2. trim the near-zero parts, see `Motion.trim` 
        3. change dt to uniform, see `Motion.change_dt`
        4. rename the motions to a series.  

        This procedure is always logged to the current folder.  

        Parameters
        ----------
        `new_name` : str  

        `run_steps` : List<int>  
            The step number to run. Users can select which one out of the 4 steps to run.
            The order will also follow this List.  

        `dt` : float  
            The uniform delta t.   

        `new_folder` : str  

        `trunc_dict` : dict  
            A dict that defines the `time` parameter in `Motion.truncate`  
            In the dict, keys are the indexes of the motions in the suite, values are the `time` values.

        `left`, `right` : float  
            The parameters of PGA to trim. Refer to `Motion.trim`.  

        Returns
        -------
        `suite` : `Suite`

        Outputs
        -------
        All the logs and the pre- and post- beautify average spectrum.
        '''
        suite = self.copy()
        spectrum_hist = suite.get_average_spectrum()
        suite.plot_all_spectrums(engine='pyplot')
        for i, motion in enumerate(self.motions):
            mo = motion.copy()
            mo.name = f'{self.name}-{i}-b'
            for step in run_steps:
                if step == 1:
                    if i in trunc_dict:
                        mo = mo.truncate(time=trunc_dict[i], new_name=mo.name+'1')
                elif step == 2:
                    mo = mo.trim(new_name=mo.name+'2', left=left, right=right)
                elif step == 3:
                    mo = mo.change_dt(dt=dt, new_name=mo.name+'3')
                elif step == 4:
                    mo.name = f'{i:02}' if len(self.motions) <= 100 else f'{i:03}'
                else:
                    print(f'step={step} is not recognized and ignored.')
            suite.motions[i] = mo
        spectrum_new = suite.get_average_spectrum()
        plt.figure()
        plt.plot(spectrum_hist.periods, spectrum_hist.spectrum, label='before')
        plt.plot(spectrum_new.periods, spectrum_new.spectrum, label='after')
        plt.plot(suite.target_spectrum.periods, suite.target_spectrum.spectrum, label='target')
        plt.legend()
        plt.suptitle(f'Beautify: {self.name} -> {new_name}')
        plt.savefig(f'{self.folder}/{self.name}-beautify.png', dpi=200)
        suite.name = new_name 
        if new_folder is not None:
            suite.set_folder(new_folder)
        return suite

    def scale(self, new_name, new_folder=None, factor_x=1, factor_y=1, scale_target_spectrum=True):
        '''Scale the motions both horizontally and vertically.  

        Parameters
        ----------
        `new_name` : str  

        `new_folder` : str, None  

        `factor_x`, `factor_y` : float  
            The factors to amplify in x axis and y axis.  

        `scale_target_spectrum` : boolean  
            Whether scale the target spectrum as well.

        Returns
        -------
        `suite` : `Suite`
        '''
        suite = self.copy()
        suite.name = new_name
        suite.set_folder(new_folder)
        for i, motion in enumerate(self.motions):
            suite.motions[i] = motion.scale(
                new_name=motion.name, factor_x=factor_x, factor_y=factor_y)
        if scale_target_spectrum:
            suite.target_spectrum = self.target_spectrum.scale(
                factor_x=factor_x, factor_y=factor_y)
        if factor_x != 1:
            suite.period_bound = [p * factor_x for p in self.period_bound]
            suite.period_samples = self.period_samples * factor_x
        suite.reset_spectrum_matrix()
        return suite

    def write_TCL(self, folder=None, filename=None):
        '''Write the suite to TCL forms for OpenSees.  
        Currently the TCL contains only one dict. 
        The keys of the dict is the motion names. Then the essential properties are stored.  
        At the same time, the accelerations for each motion is written to a text file.

        Parameters
        ----------
        `folder` : str, None  
            The folder to write the tcl files.
            If is None, use `Suite.folder`.

        `filename` :  str, None  
            The tcl file name. If is None, use `Suite.name`.

        Outputs
        -------
        A tcl file and a series of acceleration files.
        '''
        folder = self.folder if folder is None else folder
        filename = self.name if filename is None else filename.replace(
            '.tcl', '') + '.tcl'
        os.makedirs(folder, exist_ok=True)
        config = ''
        config += 'dict set motion keys {'
        names = [mo.name for mo in self.motions]
        config += ' '.join(names)
        config += '}\n\n'
        for i, mo in enumerate(self.motions):
            mo.write_accel(f'{folder}/{mo.name}', time=False)
            config += f'dict set motion {mo.name} path {folder}/{mo.name}\n'
            config += f'dict set motion {mo.name} npts {len(mo.accel)}\n'
            config += f'dict set motion {mo.name} dt {mo.dt:.4f}\n'
            config += f'dict set motion {mo.name} amp {self.factors[i]:.4f}\n\n'
        with open(f'{folder}/{filename}.tcl', 'w') as f:
            f.write(config)

    def plot_all_accels(self, folder=None, filename=None, save=True, engine='pyplot'):
        '''Plot amplified accelograms to a single file.  

        Parameters
        ----------
        `folder` : str, None  
            If is None, use `Suite.folder`.

        `filename` : str, None  
            If is None, use `Suite.name` + '_accel'

        `save` : boolean  
            If True, save the file. Else, show immediately.  

        `engine` : 'pyplot' or 'plotly'  

        Outputs
        -------
        A figure is plotted.
        '''
        folder = self.folder if folder is None else folder
        if engine == 'plotly':
            filename = f'{self.name}_accel.html' if filename is None else filename + '.html'
            fig = go.Figure()
            for i, motion in enumerate(self.motions):
                fig.add_trace(go.Scatter(
                    x=np.arange(len(motion.accels)) * motion.dt,
                    y=motion.accels,
                    mode='lines',
                    name=str(i) + motion.name,
                    text=motion.name
                ))
            if save:
                fig.write_html(f'{folder}/{filename}')
            else:
                fig.show()
        elif engine == 'pyplot':
            print('Pyplot is not supported.')
        else:
            print('ERROR: engine should be in "pyplot" and "plotly"')

    def plot_all_spectrums(self, folder=None, filename=None, save=True, engine='pyplot'):
        '''Plot all amplified spectrums in one file.  

        Parameters
        ----------
        Refer to `Suite.plot_all_accels`. 

        The filename if is None will be `Suite.name` + '_spectrum'  

        Outputs
        -------
        A figure will be generated.
        '''
        if engine == 'plotly':
            folder = self.folder if folder is None else folder
            filename = f'{self.name}_spectrum.html' if filename is None else filename + '.html'
            fig = go.Figure()
            for i, motion in enumerate(self.motions):
                fig.add_trace(go.Scatter(
                    x=motion.spectrum.periods,
                    y=motion.spectrum.spectrum * self.factors[i],
                    mode='lines',
                    name=motion.name,
                    text=motion.name
                ))
            spectrum1 = self.target_spectrum
            spectrum2 = self.get_average_spectrum()
            fig.add_trace(go.Scatter(
                x=spectrum1.periods,
                y=spectrum1.spectrum,
                mode='lines',
                name='Target',
                text='Target',
                line=dict(width=5),
            ))
            fig.add_trace(go.Scatter(
                x=spectrum2.periods,
                y=spectrum2.spectrum,
                mode='lines',
                name='Average',
                text='Average',
                line=dict(width=5),
            ))
            if save:
                fig.write_html(f'{folder}/{filename}')
            else:
                fig.show()
        elif engine == 'pyplot':
            folder = self.folder if folder is None else folder
            filename = f'{self.name}_spectrum.png' if filename is None else filename + '.png'
            plt.figure(figsize=(7, 5))
            for i, motion in enumerate(self.motions):
                plt.plot(motion.spectrum.periods, motion.spectrum.spectrum * self.factors[i], label=motion.name, c='silver')
            plt.plot(self.target_spectrum.periods, self.target_spectrum.spectrum, label='Target', linewidth=3)
            averageSpectrum = self.get_average_spectrum()
            plt.plot(averageSpectrum.periods, averageSpectrum.spectrum, label="average", linewidth=2)
            plt.legend(loc='best', fontsize='x-small')
            plt.title(f'SRSS={self.get_srss():.4f}')
            if save:
                plt.savefig(os.path.join(folder, filename), dpi=200)
            else:
                plt.show()
        else:
            print("ERROR: wrong engine name. use 'plotly' or 'pyplot'.")

    def plot_individual(self, save=True):
        '''Plot all individual ground motions and the spectrum.

        Parameters
        ----------
        Refer to `Suite.plot_all_accels`.  

        Outputs
        -------
        If save, a figure for each motion will be generated.  
        '''
        for i, mo in enumerate(self.motions):
            _, ax = plt.subplots(2, 1, figsize=(7, 5))
            ax[0].plot(np.arange(len(mo.accel)) * mo.dt,
                        mo.accel * self.factors[i])
            ax[1].plot(mo.spectrum.periods,
                        mo.spectrum.spectrum * self.factors[i])
            ax[1].plot(self.target_spectrum.periods,
                        self.target_spectrum.spectrum, linewidth=2)
            plt.suptitle(
                f'{i}: Name {mo.name}\nfactor={self.factors[i]:.4}, PGA={mo.get_PGA()*self.factors[i]:.4}, dt={mo.dt:.3f}')
            filename = f'{self.folder}/motion{i}.png'
            if save:
                plt.savefig(filename, dpi=200)
            else:
                plt.show()

    def plot_interactive(self, save=True):
        '''Plot the spectrum and the acceleration time-histories to an interactive html file.
        
        Parameters
        ----------
        `save` : boolean  
            If True, write to html file. Otherwise, show directly without saving.  
        
        Outputs
        -------
        A html file with interactive plotting sliders.
        '''
        fig = make_subplots(rows=1, cols=2)
        for i in range(len(self.motions)):
            fig.add_trace(go.Scatter(
                x=self.motions[i].spectrum.periods,
                y=self.motions[i].spectrum.spectrum * self.factors[i],
                name=self.motions[i].name,
            ), row=1, col=1)
            fig.add_trace(go.Scatter(
                x=self.motions[i].get_times(),
                y=self.motions[i].accel * self.factors[i],
                    name=self.motions[i].name,
            ), row=1, col=2)
        spectrum_average = self.get_average_spectrum()
        fig.add_trace(go.Scatter(
            x=spectrum_average.periods,
            y=spectrum_average.spectrum,
            name='average',
            line={'width': 5}
        ))
        fig.add_trace(go.Scatter(
            x=self.target_spectrum.periods,
            y=self.target_spectrum.spectrum,
            name='target',
            line={'width': 5}
        ))

        def _title(i):
            return f'Motion #{i}: factor={self.factors[i]:.2f}, duration={len(self.motions[i].accel)*self.motions[i].dt:.2f}, PGA={self.motions[i].get_PGA()*self.factors[i]:.4}'
        
        steps = []
        for i in range(len(self.motions)):
            step = dict(
                method='update',
                args=[
                    {'visible': [j//2 == i for j in range(len(self.motions)*2)] + [False, True]},
                    {'title': _title(i)}
                ],
                label=str(i),
            )
            steps.append(step)
        sliders = [dict(active=0, steps=steps, len=0.9, currentvalue={'visible': False})]
        updatemenus = [dict(
            buttons=[dict(
                label='All',
                method='update',
                args=[
                    {'visible': [True] * (len(self.motions)*2+2)},
                    {'title': f'Suite {self.name}: {len(self.motions)} motions, SRSS={self.get_srss():.4}'}
                ]
            )],
            type='buttons',
            x=1,
            y=0,
            xanchor='right',
            yanchor='top',
            pad={'t': 30},
        )]
        fig.update_layout(sliders=sliders, title=f'Suite {self.name}: {len(self.motions)} motions, SRSS={self.get_srss():.4}', updatemenus=updatemenus)
        if save:
            fig.write_html(f'{self.folder}/{self.name}.html')
        else:
            fig.show()
    
    def to_dict(self, separate_motion=True):
        '''Convert all the properties to dict.  

        Parameters
        ----------  
        `separate_motion` : boolean  
            If True, the motion details are not included into the dict.  
            User has to make sure these motions are not moved.

        Returns
        -------  
        `it` : dict  
        '''
        return dict(
            folder=self.folder,
            name=self.name,
            motions=[mo.to_dict() for mo in self.motions] if not separate_motion else [
                f'{mo.folder}/{mo.name}' for mo in self.motions],
            factors=self.factors if type(self.factors[0] == float) else [
                fact.item() for fact in self.factors],
            period_bound=self.period_bound,
            target_spectrum=self.target_spectrum.to_dict(),
            period_samples=self.period_samples.tolist(),
            srss=self.get_srss(),
            average_spectrum=self.get_average_spectrum().to_dict(),
        )

    def to_json(self, separate_motion=True):
        '''Convert the properties to json.  

        Parameters
        ----------
        Refer to `Suite.to_dict`.  

        Returns
        -------
        `json` : str
        '''
        return json.dumps(self.to_dict(separate_motion=separate_motion))

    def save(self, separate_motion=True, rewrite=False):
        '''Save the suite to a .suite file.  

        Parameters
        ----------  
        `separate_motion` : boolean  
            If True, the motion details are not included into the file.  
            User has to make sure the motions are not moved.

        `rewrite` : boolean   
            If True, all the '.motion' files will be rewritten.  
            otherwise, the motion files will be ignored.

        Outputs
        -------
        Create a .suite file to store the suite.
        '''
        with open(f'{self.folder}/{self.name}.suite', 'w') as f:
            f.write(self.to_json(separate_motion=separate_motion))
        if separate_motion:
            for mo in self.motions:
                if not os.path.exists(f'{mo.folder}/{mo.name}.motion'):
                    mo.save()
                    print(f'Motion saved to {mo.folder}/{mo.name}.motion')
                else:
                    if rewrite:
                        mo.save()
                        print(f'Motion {mo.folder}/{mo.name}.motion is written.')
        print(f'Saved to {self.folder}/{self.name}.suite')

    @staticmethod
    def read_dict(it):
        '''Read the dictionary representing the properties of the class.  
        If key 'motions' is a instance of dict, then construct motions with dict.   
        else, key 'motions' should be a list of {folder}/{name}, then construct with the .motion file.  

        Parameters
        ----------  
        `it` : dict  
            The dict generated by `Suite.to_dict`  

        Returns
        -------  
        `suite` : `Suite` 
        '''
        suite = Suite(
            folder=it['folder'],
            name=it['name'],
            target_spectrum=Spectrum.read_dict(it['target_spectrum']),
            period_bound=it['period_bound']
        )
        for i, mo_content in enumerate(it['motions']):
            if isinstance(mo_content, dict):
                suite.add_motion(Motion.read_dict(
                    mo_content), it['factors'][i])
            else:
                suite.add_motion(Motion.load(
                    f'{mo_content}.motion'), it['factors'][i])
        suite.period_samples = np.array(it['period_samples'])
        suite.target_vector = suite.target_spectrum.get(suite.period_samples)
        suite.reset_spectrum_matrix()
        return suite

    @staticmethod
    def load(filename):
        '''Load a the suite from a file.  

        Parameters
        ----------
        `filename` : str  
            File that stores the suite.

        Returns
        -------
        `suite` : `Suite`
        '''
        if not filename.endswith('.suite'):
            filename += '.suite'
        with open(filename, 'r') as f:
            return Suite.read_dict(json.loads(f.read()))

    def copy(self):
        '''Copy the instance.  
        Returns
        -------
        `suite` : `Suite`
        '''
        return Suite.read_dict(self.to_dict(separate_motion=False))
